# Imports
from coopr.pyomo import *
from coopr.opt import SolverFactory
from ReferenceModel import model
import numpy

# Solve WS for given number of sample realizations with fixed X at X_WS
numSamples = 500
numX = 5
optVal = numpy.array([0 for i in range(numSamples)])

# Choose the solver
opt = SolverFactory('gurobi')

# See the result from part c
WS_X = [9.2, 21.66, 5.00, 9.59, 5.49]

for i in range(numSamples):
    datafile = './scenariodata/Scenario' + str(i+1) + '.dat'
    instance = model.create(datafile)

    # Fix values of x at x_WS
    for j in range(numX):
        instance.X[j+1] = WS_X[j]
        instance.X[j+1].fixed = True
    instance.preprocess()

    # Solve the instance
    results = opt.solve(instance)
    print "Solve" + str(i) + "th instance"
    instance.load(results)
    optVal[i] = value(instance.TotalProfit)
    
# Point estimate
EWS = optVal[:].mean()

# Interval estimate
z_val = 1.96
EWS_var  = optVal[:].var()*numSamples/(numSamples-1)
EWS_halfwidth = z_val*sqrt(EWS_var/numSamples)
EWS_CI_lo = EWS - EWS_halfwidth
EWS_CI_hi = EWS + EWS_halfwidth
  
print EWS
print EWS_CI_lo, EWS_CI_hi

